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Abstract. We discuss results from simulations of black hole formation in failing core-collapse 
supernovae performed with the code GRID, a new open- source Eulerian spherically- symmetric 
general-relativistic hydrodynamics code. GRID includes rotation in an approximate way (1.5D), 
comes with multiple finite-temperature nuclear equations of state (EOS), and treats neutrinos in 
the post-core-bounce phase via a 3-flavor leakage scheme and a heating prescription. We chose 
the favored Kq = 220 Me V- variant of the Lattimer & Swesty (1990) EOS and present collapse 
calculations using the progenitor models of Limongi & Chief fi (2006). We show that there is 
no direct (or "prompt") black hole formation in the collapse of ordinary massive stars (8M < 
^zams ^ 1OOM ) and present first results from black hole formation simulations that include 
rotation. 
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INTRODUCTION 

Core-collapse supernova explosions are the dramatic events heralding massive star death 
in core collapse. All stars in the zero-age-main-sequence (ZAMS) mass range from 
> 8M to ~ 1OOM undergo core collapse at the end of their life, but not all core 
collapse events result in a core-collapse supernova explosion. In any star of the mass 
range under consideration, core collapse separates the stellar core into subsonically 
collapsing inner core and supersonic outer core. The collapse of the former is stabilized 
by the stiffening of the nuclear equation of state (EOS) near nuclear density. Core bounce 
occurs, launching a hydrodynamic shock into the still infalling outer core. This prompt 
shock fails to blow up the star and is forced into stall by the dissociation of accreting 
iron-group nuclei and neutrino losses from the postshock region. The shock must be 
revived for core collapse to result in an explosion. The precise mode of revival provided 
by the much sought-after supernova mechanism is still uncertain (e.g., [1, 2, 3, 4, 5, 6]), 
and in a finite (though unknown) fraction of massive stars, it must fail completely to 
drive an explosion, or, at least fail to unbind the entire star so that significant fallback 
accretion occurs [7]. In the former case, a black hole (BH) inevitably forms within a 
few seconds as the stellar mantle accretes onto the protoneutron star (PNS; Collapsar 
type I [8]), while in the latter, a neutron star initially survives but may be pushed over 
its mass limit by fallback accretion (Collapsar type II [8, 9]). Both are considered as 
potential scenarios leading to a long gamma-ray burst (GRB) [8, 10]. 

BH formation in the core collapse context has been studied with full general- 
relativistic (GR) Boltzmann neutrino radiation-hydrodynamics codes in spherical 



symmetry without rotation (e.g., [11, 12, 13] and references therein), but due to the 
computational cost of these high-fidelity simulations, only few models have been sim- 
ulated. Multi-dimensional simulations of BH formation remain to be performed with 
appropriate microphysics and progenitor models, but see [14] for an exploratory work 
with simplified physics. 

We have recently implemented the new code GRID, an Eulerian spherically- 
symmetric code for stellar collapse and BH formation [15]. GRID is open-source and 
may be downloaded from http: / /www. stellarcollapse.org. GRID, in its 
present version, does not implement the full radiation transport formalism, but instead 
relies on an approximate, yet extremely computationally efficient leakage scheme 
for neutrino transport (e.g., [16, 17]) and a simple prescription for neutrino heating. 
Moreover, GRID implements an approximate way of including rotation in spherical 
symmetry a variant of which is used in stellar evolutionary calculations (e.g., [18]). 
GRID'S computational efficiency allows us to perform hundreds of model simulations 
in little time to explore systematically the conditions for BH formation in massive stars 
and the parameter space in ZAMS mass, metallicity, and rotation where BH formation 
may be the dominant outcome of stellar collapse. Detailed results of such an extensive 
study will be reported in [19]. 

In the following, we discuss some details of the GRID code, then present simulation 
results highlighting the fact that any BH forming core collapse passes through a PNS 
phase and is never "direct". We then show first results from simulations of rotating BH 
formation. 



GRID follows the 3 + 1 approach to numerical relativity, slicing 4D spacetime into 3D 
spacelike hypersurfaces along a timelike normal (e.g., [20]). This introduces two gauge 
quantities, the lapse function a and the shift vector j8 whose choices are not a priori 
fixed. The lapse controls how time changes between two consecutive slices while the 
shift describes how coordinates change from one slice to the next. In GRID, we adopt 
spherical symmetry and the polar-slicing, radial-gauge choice (e.g., [21]) on an Eulerian 
grid, resulting in j8 = and a Schwarzschild-like invariant line element of the form 
(assuming c — G — M = 1 here and in the following), 



where the lapse a and the metric function X are functions of a metric potential 4>(r,f) 
and of the enclosed gravitational mass M grav (r,f) = m{r,t\ 



THE GRID CODE 



ds 2 = -a(r,t) 2 dt 2 + X(r,t) 2 dr 2 + r 2 d£l 2 , 



(1) 





(2) 



The detailed form of 4>(r,f) can be found in [15, 21]. Here we point out only that the 
choice of <E> ensures that a is singularity-avoiding and, hence, drops to very small values 
when a physical singularity forms, thus minimizing its evolution with coordinate time. 



It is important to note that the metric function X becomes singular at r — 2M and thus 
forbids numerical evolution beyond BH formation in this gauge (a feature in common 
with many other ID codes [11, 13]). 

GRID'S hydrodynamics module follows the flux-conservative Valencia formulation 
of GR hydrodynamics in the form of [22]. The scheme is semi-discrete in space and is 
finite- volume with piecewise-parabolic reconstruction of interface values and employs 
the HLLE Riemann solver [23]. Time discretization is handled via 2nd- or 3rd-order 
Runge-Kutta integrators using the Method of Lines [24]. 

A specialty of GRID is its approximate inclusion of rotation in spherical symme- 
try (1.5D). This is accomplished by solving an advection equation for the angular mo- 
mentum and including an angularly averaged centrifugal term in the radial momentum 
equation, in the Lorentz factor, and in curvature terms to account for centrifugal force, 
angular momentum flux, and rotational energy [15]. 

GRID operates with a general EOS interface and has been tested with variants of the 
Lattimer-Swesty EOS [25] and the H. Shen EOS [26] tables of which we make available 
for download at http : //stellarcollapse . org and describe in [15]. 

Deleptonization and neutrino transport are handled by GRID differently in the pre- 
bounce and postbounce phases. In the former, we employ the simple parameterization of 
the electron fraction Y e as a function of rest-mass density p put forth by [27]. In the lat- 
ter, we use a 3-flavor (v e , v e , and v x = { v M , , v T , v T }) energy-averaged leakage scheme 
constructed along the lines of [16, 17] and described in [15]. It yields energy-integrated 
neutrino luminosities that are within ~ 20% of those predicted by full Boltzmann trans- 
port calculations. Neutrino heating by charged-current absorption of v e and v e is handled 
via 



where L Vi (r) is the luminosity as set by the energy leakage rate interior to radius r, 
Cheat, Vi i s the energy- averaged absorption cross section, m u is the atomic mass unit, X; is 
the mass fraction of the absorbing particle (proton or neutron), and (1/F y .) is the mean 
inverse flux factor which we approximate analytically as a function of the optical depth 
T by comparing to angle-dependent radiation transport calculations [3, 15]. The factor 
e~ 2Zy i is used to attenuate heating at high optical depths where neutrinos and matter are 
in equilibrium, /heat is a scaling factor that may be used to dial in higher heating rates, 
but is set to 1 here, causing explosions to fail in all models. Once the heating rate is 
computed, the luminosity is reduced accordingly to maintain energy conservation. 



NO DIRECT BLACK HOLE FORMATION IN CORE COLLAPSE 



The idea prevails that stars in the upper half of the ~ 8M to ~ 1OOM mass range 
collapse to BHs directly. This is incorrect, but, by an unfortunate choice of words in 
recent work [8, 28] that can easily be misunderstood, has become widely believed. 

According to Thome's Hoop Conjecture [29], a BH forms when a given amount of 
mass-energy collapses through its own Schwarzschild radius R$ = 2M (in geometric 
units), hence must achieve a compactness M/R > 0.5. If prompt BH formation were 
to occur in core collapse, the collapsing core (or parts of it) must reach this level of 
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FIGURE 1. Metric coefficient X = 1 / a/1 — 2m/r as a function of enclosed baryonic mass at the time 
of core bounce and greatest inner-core compactness. We plot X for GRID collapse simulations of models 
from [30] with ZAMS masses from 20 to 12OM . If a BH were to form, X would diverge. Due to the 
universality of core collapse, inner core masses vary only little with progenitor and the compactness at 
bounce stays moderate independent of progenitor mass / model. There is no prompt formation of BHs. 
While GRID'S treatment of deleptonization during collapse is approximate, it agrees well qualitatively 
and quantitatively with more quantitatively accurate simulations (e.g., [27]). 

compactness prior to or at bounce. In the following, we demonstrate that this is never 
the case for stars in the above mass range and explain why. 

Goldreich & Weber [31] and Yahil [32] demonstrated analytically what has been 
confirmed numerically numerous times: In collapse, the iron core separates into the 
homologously (v oc r, in sonic contact) infalling inner core (IC) and the supersonically 
collapsing outer core. It is only the inner core that plunges to nuclear density and 
significant compactness in the final phase of core collapse and it is essentially its mass 
(Mic) that must be pushed below its Schwarzschild radius to make a BH. 

Mic at bounce is proportional to Y^ lc (l + r/sfc) [33], where r/ ^ 0.1 and where Y e ^c 
and sic are the mean Y e and specific entropy of the inner core at bounce, respectively. 
Both quantities are coupled through electron capture, neutrino transport, and the equal- 
isation of neutrinos and matter above trapping density. The physics governing Y e jc 
and sic during collapse is general and independent of the conditions prevailing in a 
particular collapsing star [27, 1]. Hence, core collapse is universal, rather independent 
of initial conditions, and Mic falls into the range of ~ 0.4 — O.6M in the nonrotating 
case 1 [35, 34]. The Hoop Conjecture would require this mass to be compressed into 
1.3 — 1.8km, which even for the softest plausible nuclear EOS (e.g., the Kq = 180 MeV 
variant of the Lattimer-Swesty EOS [25, 36]) does not occur before or at bounce. 

In Fig. 1, we plot the metric coefficient X (Eq. 2) at the time of core bounce and as a 
function of enclosed baryonic mass. X is shown for a set of solar-metallicity models 
in the 20 — 12OM mass range. These are drawn from the Limongi & Chieffi [30] 



1 Rotation can increase the size of the homologous region and, hence, Mic significantly, but also limits 
inner core compactness [34]. 
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FIGURE 2. Example results from GRID simulations of a 4O-M progenitor of [30] set up to spin 
with asymptotic specific angular momenta ji$ = {0,1,2} (see text for details). Bottom: Evolution of 
the central density as a function of postbounce time. Rotational support leads to a lower accretion rate, 
slower PNS contraction, and a larger maximum PNS mass. At BH formation, the PNS baryonic mass 
is 2.37, 2.39, and 2.46M , for the ji6 = 0, 1, and 2 cases, respectively. Top: Ratio of rotational kinetic 
energy T to gravitational binding energy \W\ in the two spinning simulations. In both models, a local 
maximum is reached at bounce and is followed by a secular increase in the postbounce phase as the PNS 
core contracts and spins up. The = 2 run reaches values of r/|W| that may make it susceptible to a 
secular nonaxisymmetric instability [38] before BH formation. Both models may develop a low -T/\W\ 
corotation-type nonaxisymmetric instability (e.g., [39]). 



model set (using the mass-loss rates of [37]) and collapsed with GRID without rotation, 
using the Lattimer-Swesty K = 220 MeV (LS220) EOS, and an analytic fit to the 
deleptonization trajectory Y e (p) obtained by [27]. The figure clearly underpins the 
universality of core collapse and demonstrates that none of the collapse models come 
even close to BH formation before or at bounce (for which X 3> 1). Variations of 
the X(m) profiles with progenitor mass are due primarily to larger inner core masses 
owing to higher temperatures/entropies in more massive progenitors. In a full radiation- 
hydrodynamics treatment, one would expect these variations to be even smaller due to 
increased deleptonization in hotter cores. 



We perform 1.5D rotating collapse simulations of the 4O-M progenitor of [30] with 
the LS220 EOS. Rotation is added when the peak collapse velocity of the core reaches 
1000 kms -1 via the rotation law Q(r) = Qq(1 + r 2 /A 2 ) -1 . A is set to the radius at which 
the enclosed baryonic mass is 1M , resulting in nearly uniform rotation throughout 



THE FORMATION OF SPINNING BLACK HOLES 



the inner core [40]. We perform three calculations with = {0,0.9, 1.8}rads , cor- 
responding to specific angular momenta at infinity in units of 10 16 cm 2 s _1 of 715 «, = 
{0,1,2}. 

In the lower panel of Fig. 2, we show the time evolution of the central density in 
the three calculations as a function of postbounce time. The nonrotating model forms 
a BH at ~ 1.14 s and has a maximum baryonic (gravitational) PNS mass of ~ 231 M® 
(rsj 2.19M ). In the rotating models, centrifugal support reduces the accretion rate and 
slows down the contraction of the PNS. This alters the PNS structure and leads to an 
onset of PNS collapse at later times (scaling roughly with Qfy and at lower central 
densities. The times of BH formation are ~ 1.19 s and ~ 1.36 s in the 7*15,00 = 1 and 
7x6,00 = 2 model, respectively. Interestingly, rotational support has only a small effect on 
the maximum PNS masses in the models presented here. This is due primarily to the 
fact that the PNS cores in our models are uniformly spinning. A significant increase of 
the maximum PNS mass is expected only in differentially rotating PNSs [41]. We find 
baryonic (gravitational) PNS masses at BH formation of ~ 2.39M (~ 2.22M ) and 
~ 2.46 M® (~ 2.3OM ) in the two spinning models. 

The top panel of Fig. 2 depicts the time evolution of the ratio of rotational kinetic 
energy T to gravitational energy \W\ in the spinning models. In both, r/|W| increases 
secularly after bounce as the PNS accretes and contracts and reaches a maximum value 
before the onset of PNS collapse of ~ 0.20 (~ 0.05) in the 716,00 = 2 (716,00 = 1) model. 
Both models stay below the threshold r/|W|d yn ~ 0.27 for the high-r/|W| dynamical 
nonaxisymmetric instability (e.g., [42]). The rapidly spinning PNS of the 716,00 = 2 
model stays above r/|W| sec ~ 0.14 for ~ 1 s, a time likely sufficiently long for a secular 
(gravitational- wave or viscosity driven) nonaxisymmetric instability to arise in 3D (e.g., 
[38]), redistributing/radiating angular momentum and thus effectively limiting the PNS 
core spin. In addition, both models may be susceptible to shear instabilities that may 
generate nonaxisymmetric structure (e.g., [39]) and/or magnetic flux (e.g., [43]) and 
redistribute angular momentum. In our present 1.5D calculations such processes and 
instabilities are absent and the BH formed in the rapidly spinning 716,00 = 2 model has 
an initial spin parameter a* — J/M^ mw of ~ 0.81. Its more slowly spinning 716,00 = 1 
counterpart forms a BH with a* ~ 0.41. 

DISCUSSION 

In this contribution to the proceedings of the 10th International Symposium on the 
Origin of Matter and Evolution of the Galaxies (OMEG10), we have highlighted results 
of a small set of stellar collapse and BH formation simulations carried out with the 
new open-source 1.5D code GRID. We numerically demonstrated what has long been 
known on the basis of analytic arguments [31, 32], namely that a BH is never formed 
promptly in the core collapse of ordinary massive stars with masses between ~ 8M and 
~ 1OOM . Any core collapse event, if ultimately resulting in BH formation or not, passes 
through a protoneutron star phase in which neutrinos and, quite likely, gravitational 
waves, are emitted for at least a few hundred milliseconds. Direct subsidence into a BH is 
possible only in much more massive stars that become radially unstable before forming 
a hydrostatic iron core (e.g., [44]), provided they do not experience a pair-instability 



driven thermonuclear disruption instead of collapse. 

Massive stars, in particular those favored as the progenitors of long GRBs, may have 
rapidly spinning cores [10] that will make spinning BHs. Using GRID'S 1.5D rotation 
feature, we have performed a first small set of spinning BH formation simulations. Rota- 
tional effects increase the time to BH formation significantly (and, thus, allow a potential 
explosion mechanism more time to operate!), but, due to almost uniform rotation in the 
PNS core, have a much smaller impact on the maximum PNS mass. We also find that the 
rotation rate r/|W| of PNSs increases significantly during the postbounce accretion and 
contraction phase. PNSs with early postbounce values of T/\W\ below the thresholds 
for secular or dynamical rotational instabilities may surpass these within a few hundred 
milliseconds. Nonaxisymmetric deformation and the associated angular momentum re- 
distribution and/or emission of angular momentum in gravitational waves may put a 
natural limit on the maximum PNS spin and, in consequence, on the birth spin of BHs. 

Stellar-mass BH formation is a process that occurs frequently in the universe and 
core collapse is its natural site. As observations suggest ([? ] and references therein), 
it may be the generic ultimate outcome of core collapse in stars more massive than 
~ 20 M . In the context of the core-collapse supernova - GRB connection, BH formation 
is a necessary ingredient for the collapsar scenario to work (e.g., [8, 10]). The current 
theoretical understanding of core-collapse supernovae and BH formation is still only 
partial. Ultimately, it will be necessary to establish a firm quantitative mapping between 
ZAMS conditions (mass, metallicity, angular momentum) and the outcome of stellar 
collapse. This is a problem as much in stellar evolution as in core-collapse supernova 
theory and will require advances in both fields. 

With the GRID code we have created an open-source tool that allows us and others to 
study BH formation in failing core-collapse supernovae and to investigate its systematics 
and the characteristics of nascent BHs with variations in presupernova stellar structure 
and rotational configuration. While a quantitatively robust mapping to ZAMS conditions 
may not be possible on the basis of currently available stellar evolutionary models, robust 
qualitative features and trends can be derived in parameter studies [19]. 
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